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I. INTRODUCTION 

Density functional theory (DFT) in general and the 
Kohn-Sham scheme in particular are among the most im- 
portant tools for electronic structure calculations. This 
success is based on a rigorous foundation in terms of 
the Hohenberg-Kohn theorem on the one handi, and in- 
creasingly sophisticated approximations to the exchange- 
correlation functional E xc on the other. The local den- 
sity approximation (LDA)2 proved to be of far greater 
applicability than originally expected. Its great advan- 
tages are simplicity and reliability in the sense that 
its shortcomings are qualitatively predictable, but the 
price to be paid is limited accuracy. Generalized gra- 
dient approximations (GGA's)^i^ improved accuracy 
to a level that made DFT useful for many chemical ap- 
plications. But despite their many successes, GGA's 
cannot be regarded as the final stage of functional 
developments^. Further improvements in accuracy are 
expected from functionals that include partial or full ex- 
act exchange. Globa&2ii2iiiii2 and locali^ hybrids and 
hyper-GGA'sii^ fall into this class. 

Using the exact exchange energy in DFT calculations is 
promising from many points of view. Full exact exchange 
cancels the spurious Hartree self-interaction energy, cur- 
ing in a systematic way one of the most notorious DFT 
problems. It leads to the correct high-density limilji^, 
which is exchange dominated. Great improvements in 
the Kohn-Sham eigenvalue spectrumi§*i£, semiconductor 
bandstructure and excitationsi&i2i22i2I and nonlinear op- 
tical properties 2 ^ have been reported using the exact ex- 
change energy. The exact exchange functional also leads 
to the correct asymptotic (r — > oo) behavior of the Kohn- 
Sham potential, which turns out to show surprising fea- 
tures not present in any of the common approximations, 
as discussed in Refs. Il6l23l and below. Finally, the idea 



of incorporating one more known exact ingredient into 
the energy functional is appealing in its own right. 

But the goal of using the exact exchange energy or, 
more generally, any orbital functional sclf-consistently 
in Kohn-Sham calculations raises the question of how 
to construct the corresponding exchange-correlation po- 
tential. E xc is still implicitly a functional of the den- 
sity, but explicitly only its dependence on the set of 
occupied Kohn-Sham orbitals is known. Thus, calcu- 
lating the corresponding exchange-correlation potential 
^xc( r ) = 3E xc [n]/6n(r) cannot be done straightforwardly 
by explicitly taking the functional derivative with respect 
to the density n(r). Instead, the potential must be cal- 
culated from the optimized effective potential (OEP) in- 
tegral equationSii2Si2&. This approach is size-consistent 
if the orbital functional is invariant under unitary trans- 
formations of the occupied orbitals and is not too radi- 
cally nonlocalSi. Size consistency is achieved for exact ex- 
change, hybrid functionala 8 i 9 i 10 i n i 12 i 13 , meta-GGA's^a 
and hyper-GGA'sLii. 

Due to its complexity, direct solutions of the 
OEP integral equation so far have been restricted 
to effectively one-dimensional systems with spherical 
symmetr y 2 ^??:? 1 !? 2 :??!? 4 !?^?^ , and the calculations for fi- 
nite systems have been based on the program developed 
in Ref. |25[ Hopes are low to generalize this approach 
to higher dimensions. The OEP for three-dimensional 
systems has been constructed by directly evaluating the 
response function using basis se tsi2ii2i2Li22iSi. This ap- 
proach has proven successful, but requires considerable 
technical expertise: Summing over not only occupied but 
also unoccupied Kohn-Sham orbitals is required, and the 
necessary inversion of the response function can be cum- 
bersome. And while the OEP total energy can be ac- 
curately obtained, the potentials and Kohn-Sham eigen- 
values for finite systems may suffer from basis-set lim- 
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itations. A detailed discussion of the m ethod and its 
limitations can be found in Refs. l4d4ll42L Alternatively, 
the short-range part of the potential itself has been ex- 
panded in a basis set and the expansion coefficients cho- 
sen to minimize the total energy^^^. This approach is 
appealing because of its directness. But the — e 2 /r-decay 
that the OEP is thereby forced to take everywhere in the 
asymptotic region of a finite system has to be reconsid- 
ered in the light of new conclusions about the asymptotic 
behavior of w x (r) presented in Refs. Il6ll23l and discussed 
below. Also, the method is hard to implement within 
fully numerical schemes for electronic structure calcula- 
tions, and those continue to grow in importance^SiiLiS. 

In this manuscript we discuss how the OEP can eas- 
ily be constructed from the first-order orbital shifts that 
have been introduced into OEP theory to justify the KLI 
approximation^. We first present in Section [H] a new 
proof for the OEP equation that explains why this com- 
plicated equation can be cast in a simple form. In Section 
IIIII we compare and contrast Hartree-Fock theory and 
exact-exchange OEP. Different methods to construct the 
OEP from the orbital shifts are analyzed in Section ITvl 
together with indicators for the accuracy of an OEP so- 
lution. The exact-exchange OEP for spherical atoms and 
three-dimensional sodium clusters is calculated. Its long- 
range asymptotic behavior is investigated in Section Ivl 



We show that the behavior lim r ^ 00 w x (r) 



7r + C, 



with C ^ on nodal surfaces of the highest occupied 
Kohn-Sham orbital, is a common situation for metal clus- 
ters. Finally, in Section lVH we calculate the static electric 
polarizability for the neutral, even-electron clusters Na2 
- Naio in the exchange-only approximation, and compare 
LDA, KLI and exact OEP results. 



is the one that makes the change in the density vanish to 
first order in the perturbation 



A%(r) = u xci(T (r) - u xco .(r), 



(3) 



when this perturbation is applied to the Kohn-Sham sys- 
tem. The ^icr(r) in Eq. are the Kohn-Sham orbitals, 
i.e., the solutions of 



(h KSa - Eia^j tp ia (r) = 0, 



(4) 



where fiKSa- = — (7i 2 /2m)V 2 + i>KSer( r ) is the Kohn-Sham 
Hamiltonian. The self-consistent Kohn-Sham potential 



^KSa(r) = u(r) + un(r) + v XC(7 (r) 



(5) 



is the sum of the external potential v(r), the Kartree 
potential t>H( r ) = / d 3 r'e 2 n(r')/\r — r'\, and the 
spin-dependent exchange-correlation potential u XCCT (r) = 
SE xc [n]/Sn a (r). (v(r) can also depend upon a, but usu- 
ally does not.) ipia(r) is the negative of the first-order 
perturbation-theory shift that results if <fi a {^) is sub- 
jected to the perturbation of Eq. |JS}, i.e, 



-V4( r ) 



(0) 



y J gUP [w(rQ - iw(r')] ^(r') dV . 

For the sake of simplicity, we here assume non-degenerate 
orbitals. The extension to the degenerate case is unprob- 
lematic: As shown in Appendix B of Ref. the restric- 
tion j ' ^ i in the sum is to be replaced by Sj a ^ e^. 
Note that, if the exact exchange energy 



II. THE OEP EQUATION IN TERMS OF THE 
FIRST-ORDER DENSITY SHIFT 



As shown by Krieger, Li and Iafrate^S and further clar- 
ified by Grabo, Kreibich, Kurth and Gross^, the OEP in- 
tegral equation for the spin-dependent exact Kohn-Sham 
exchange-correlation potential v xca (r) can be written in a 
form that takes a very simple interpretation. At the end 
of a long argument involving repeated application of the 
chain rule for functional derivatives and linear response 
theory, the OEP equation is written in the form 



(1) 



Eq. Q says that the optimum (i.e., yielding the low- 
est Kohn-Sham energy) potential v xca (r) to replace the 
orbital-dependent potential 



^xczer (f) 



1 6E XC [M] 



(2) 



[M] 



■E 



^L( r M CT ( r >ia(r)^ CT (r') 
— : d r d r. 



r — r' 



(7) 

is substituted for E xc , then Eq. J2J j us t yields the orbital- 
dependent Kartree-Fock potential^ 



r(r) = - 



Eq. Q is intuitively appealing for the following rea- 
sons: (1) It is a statement about the electron spin den- 
sity, the central quantity of density functional theory. 
Going over from the optimized effective potential to the 
orbital-dependent potentials should not change the den- 
sity much. In important limits (one- and two-electron 
ground states, and the uniform electron gas), the density 
should not change at all. (2) Once the OEP orbitals and 
orbital energies are fixed, Eq. is homogeneous of de- 
gree one in the perturbation of Eq. , making v xca (r) a 
kind of average of the M xf yr(r). A simple average was pro- 
posed in Eq. (75) of Ref. |5(1 and averaged potentials are 



3 



popular in time-dependent DFT— . A sophisticated but 
still approximate average was proposed in Ref. I30L (3) 
Eq. contains just enough information to define v xca (r) 
uniquely to within an additive constant, given the OEP 
orbitals and orbital energy differences: Imagine discretiz- 
ing the r space into a collection of M points ri, ...,tm, 
and solving Eq. JJJ as a set of M — 1 linear equations for 
the M — 1 unknowns v xca (ri), ...,v xca (rM-i)- Because 
the number of electrons cannot change, only M—l equa- 
tions are linearly independent, leaving v xc „(ym) as a free 
parameter. This parameter is normally chosen so that 
]im I ._ >0O u XC(r (r) = 0. 

The OEP equation has previously been investigated 
from different perspectivea 24 i 2; ?'??i?^?Ti 4 ?i i P 2 i ; P?ii? 4 . But an 
equation as charmingly simple as Eq. Q also deserves 
a simple proof. We provide one in the following. For 
notational simplicity, we drop the spin index a in the 
rest of this section. 

We start from a given expression for the total energy 
in terms of some set of orbitals {(/>}, 

S?K#] = 2^ E / \ V Mr)\ 2 d 3 r + J V (r)n(r)d 3 r 



i=l 



e 

T 



n W n ^ r > d 3 r + aExc[m (9) 



r — r 



where 



as those normalized linearly-independent orbitals that 
yield a given density via Eq. H1(J|) and deliver the low- 
est extremum of E"[{4>}}. (This definition is discussed 
in Appendix The density n orb ' Q (r) that minimizes 
E° Yh ' a \n] results from the orbitals that are solutions of 



2 m 



V 2 + W (r)+« H (r)+< CJ *(r)Ur b ' Q (r) 



orb, a orb, a/ \ / A A \ 

e< V l (r). (14) 



The requirement of a local, orbital-independent potential 
in the definition of E® EP ' a {n] can be understood as an 
additional constraint not present in E° lb,a [n). Therefore, 
if we define a functional A a [n] by 



E" 



E? EP ' a \n] 



A c 



(15) 



we know that A a [n] < 0. From the minimum principle for 
E° rh ' a [n), and from the preceding statement, it is clear 
that 



E° 



"] < E° Ib - a [ 



b,ar OEP, a 



}<E 



OEP.ar OEP,ai 



(16) 

For a given orbital functional, one achieves a lower energy 
with the orbital-dependent potentials than with the OEP. 

Eq. |JJ| can now be deduced directly from the a- 
dependence of the two energy functionals. We first ob- 
serve that 



n 



n(r)=J2Mr)\ 



(10) 



„OEP,q 



(N;r) 



Sn(r) 



(17) 



The real parameter a is introduced for the sake of later 
arguments and can be set to 1 at the end. Based on Eq. 
© , two different density functionals can be defined. The 
Kohn-Sham (i.e, OEP) functional is defined by choosing 
the {(/>} to be the Kohn-Sham orbitals {ip}. This means 
that for a given density n(r), we find a local potential 
v s (r) such that the solutions of the independent-particle 
Schrodinger equation 



2m 



+ v s {r) ipi(r) = enpi(r) 



(11) 



yield the given density via Eq. I|10|l . The orbitals are 
clearly functionals of v s (r) which is itself a functional of 
n(r). The Kohn-Sham energy functional is then 



e: 



OEP.a 



"]=»]}]■ 



(12) 



Eq. I|12|) is minimized by the density n > a (r), and this 
density is obtained by choosing for v s (r) the potential of 
Eq. JSJ) (with u xc multiplied by a, see below). We define 
a second, "orbital" density functional 



orb, a 



n}}} 



(13) 



as the energy obtained by evaluating Eq. JHJ with a dif- 



ferent set of orbitals {ip mh ' a }. The (p° rh ' a are defined 



and 



<«:(•!>}; r ) = a 



1 SE XC [{0}] 



(18) 



i.e., both the orbital- independent and the orbital- 
dependent potential depend linearly on a. We further 
observe that, in contrast to the ip° lh ' a (r), the Kohn-Sham 
orbitals for a given density do not depend on a: The 
Hohenberg-Kohn theorem tells us that, for a given den- 
sity, v s (r) and thus the orbitals <Pi(r) are uniquely fixed. 
Finally, we observe that for a — 0, the functionals Eq. 
(|T^|l and Eq. l(T3|l are trivially identical. (Thus, in Eq. 
G3 below, E% lh > [n] = E° EP -°[n}.) 

From the second observation and Eq. @ one deduces 
that, for fixed external potential v(r) and fixed density 
n(r), the a-dependence of the OEP-functional is 

£° Ep .>] = E™ p >°[n] + aE xc [{<p[n]}}. (19) 

In contrast, -E° rb,Q [n] in general will not have such a sim- 
ple dependence on a since the (p° lh ' a for a given density 
will change with a. Thus, the power series 



E^[n] = E^[n]+aE^[n] 



a 2 E° 



(20) 



where formally E°^[n] = (dE° lh ' a [n]/da)\ a=Q , etc., will 
have non-vanishing terms beyond the one linear in a. 



4 



However, the crucial fact to note is that the term linear 
in a is the same in E° EP ' a {n] and E° rh ' a [n}: At a = 
(Hartree approximation), the orbitals are the same, i.e., 
ip° v ,0 (r) = <Pi(r), and we know that there can be no first- 
order change in E° rh,a [n] due to a change in the orbitals, 
because we defined the <p° T '"(r) as extremizing orbitals. 
Therefore, E°$[n] = E xc [{<p[n]}}. 

This means that the functional A a [n] of Eq. I|15|l van- 
ishes to linear order in a. It is then intuitively clear that 
the densities that minimize E° rh ' a [n] and £'^ )EP,Q! [n] must 
also be the same to linear order in a. In other words, 
the difference to first order in a between n orb,a (r) and 
n P ' Q (r) must vanish, and, because of Eqs. H17|) and 
(|f 8fl . this difference is just given by the left-hand side of 
Eq. Q). A detailed version of these arguments is pre- 
sented in Appendix 151 

Note that the determinant of orbitals satisfying Eqs. 
(|II|I or (|14|) and achieving the lowest total energy need 
not necessarily satisfy the Aufbau principle, but might 
have unoccupied orbitals with orbital energies lying be- 
low occupied ones. Note further that, in the absence of a 
magnetic field, all the orbitals can be chosen to be real. 



III. HARTREE-FOCK VS. EXCHANGE-ONLY 
OEP 

As an important example of the ideas of Section [HI 
let E xc [{(/>}} in Eq. ^ be the Fock exchange integral of 
the occupied orbitals of Eq. (0). Then the lowest ex- 
tremum over linearly-independent orbitals in the defini- 
tion of E° lh ' a [n] becomes the minimum over orthonormal 
orbitals. For given external potential v(r) and electron 
number N, the </?° rb ' Q_1 (r) are the canonical Hartree- 
Fock orbitals. £° rb < Q=1 [ n] is the density functional for 
the Hartree-Fock energy defined by the Levy constrained 
searchS^ (but not known as an explicitly calculable 
functional of the Kohn-Sham orbitals or of the density). 
Its minimizing density n orb ' Q=1 (r) is the Hartree-Fock 
density, and E° lh ' a=1 [n orb < Q=1 ] is the Hartree-Fock en- 
ergy. The minimizing problem for this functional can be 
formulated within Kohn-Sham theory, but the resulting 
Kohn-Sham orbitals (pi([n° rh ' a=1 ]; r) are not the Hartree- 
Fock orbitals, although both sets of orbitals yield the 
same density. The relation between Hartree-Fock theory 
and DFT has also been studied in Ref. |5?1 

On the other hand, E® EP ' a [n] is by definition^ the ex- 
act exchange-only density functional of Kohn-Sham the- 
ory (an explicitly calculable functional of the Kohn-Sham 
orbitals, and one which defines the lower limit of the in- 
tegrand in the coupling constant or adiabatic connection 



formula for E x 



158.59 



). The inequality following Eq. 



(|I5|I shows that the Hartree-Fock exchange energy for a 
given density includes a negative contribution which from 
the exchange-only OEP perspective is a small part of the 
correlation energySSiSi. (The local density approxima- 
tion to this part of the correlation energy vanishes, while 
its second-order gradient coefficient diverges 2 ^, making it 



very hard to approximate.) Since e 2 plays the same mul- 
tiplicative role as a, the exact exchange energy of Kohn- 
Sham theory2& for a given density is purely of order e 2 , 
while the Hartree-Fock exchange energy for the same den- 
sity agrees with it to order e 2 , differing in order^S e 4 . 
Eq. Ultjfl also shows that the Hartree-Fock total energy 
must be less than the total energy of the exact exchange- 
only OEP theory. These exact-exchange conclusions are 
not new2&, yet these subtle distinctions continue to cause 
some confusion. 

The Hartree-Fock and exchange-only OEP minimizing 
densities for a given v(r) also agree to order e 2 , differing 
from one another and from the correlated ground-state 
density in order e 4 . Since the highest-occupied orbital 
energy controls the asymptotic decay of the density, the 
highest occupied orbital energies must also agree to or- 
der e 2 with one another, and with minus the correlated 
first ionization energySSiS^. This helps to explain the 
observed^ but previously unexplained (except by Eq. 
(|24|l below) closeness of these quantities. 

Given v(r), the Hartree-Fock and exchange-only OEP 
orbitals and other orbital energies differ in order e 2 . A 
good approximation to the occupied Hartree-Fock or- 
bitals can be achieved by a unitary transformation of 
the occupied canonical OEP orbitals 41 . In particular, 
after mixing in the highest occupied OEP orbital, all 
the canonical Hartree-Fock orbitals so approximated will 
decay asymptotically with the same exponent. How- 
ever, exact agreement is not expected even to order e 2 , 
since Eq. 10 shows that to order e 2 each exact canon- 
ical Hartree-Fock orbital is a superposition of all the 
exact-exchange occupied and unoccupied OEP orbitals. 
Thus, for a given v(r), the Hartree-Fock and exchange- 
only OEP Slater determinants (which are invariant un- 
der unitary transformations of their occupied orbitals) 
also differ to order e 2 , and only the Hartree-Fock de- 
terminant agrees to order e 2 with the singles-only con- 
figuration interaction expansion of the correlated many- 
electron wavefunction (Brillouin's theorem^). 



IV. ITERATIVE SOLUTION OF THE OEP 
EQUATION 

A. Constructing the OEP from the orbital shifts 

The fact that the OEP can be expressed solely in terms 
of the occupied Kohn-Sham orbitals and their first-order 
shifts leads to a simple scheme for its construction. As a 
straightforward result from first-order perturbation the- 
ory, the partial differential equations 

(h K Sa ~ £i CT )V4( r ) = 

-bxctr(r) - u xcia {r) - (v xci(7 - u XC i a )](p* a {v) (21) 

for the orbital shifts ^>i CT (r) are obtained, where 
v X cia = J <^* CT (r)u XCCT (r)v? io .(r)d 3 r and u Kcia = 
J v5* .(r)M XC i IT (r)(/j i(T (r) d 3 r. (Note that the same equa- 
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tions can be obtained by other arguments'^). By solv- 
ing Eq. (|2*T1) for v^sa^*,?, inserting the result into Eq. 
JU multiplied by UKSo-( r ), and solving for v xca , the exact 
exchange-correlation potential can be expressed explic- 
itly in terms of only the occupied Kohn-Sham orbitals 
and their shift o 30 ! 36 : 



*w(r) = 

1 — 1 



2m 



+ c.c. 



(22) 



Invoking Eq. (JJJ once more, the last term can be rewrit- 
ten to obtain the expression from which the KLI approx- 
imation - neglecting the terms involving the (r) - has 
been justified as a mean- field approximation 30-36 : 



73 H {W( r )| 2 Kc«r(r) + ( 



2n CT (r) ^ 



^xcicr ^xcicr J 



-V-[V* ff (r)V^(r)] +c.c 



(23) 



Since the potential v xca (r) of Eq. I|23(l is fixed only up to 
an additive constant, one is at liberty to choose one of 
the constants v XC i a — u xcilJ freely. The usual choice i a 30 i 36 



0. 



(24) 



to make the potential vanish at infinity. We will discuss 
this point in greater detail in Section Ivl 

Together with the orthogonality condition 
j ipf a (r)(pi a (r) d 3 r = that follows from Eq. JBJ, 
Eq. 121fl uniquely determines the orbital shifts. In 
case of degeneracies, the orthogonality condition reads 
/ ^* (T (r) l £ifl(r) d 3 r = for all pairs with e ia = e Jcr . 
In Ref. |23 it was demonstrated that the ■0i* (T ( r ) can 
easily be calculated from Eq. 121|) and used to construct 
the OEP. The basic idea is to combine the Kohn-Sham 
equations Q with the orbital shift equations (|21|) in a 
self-consistent iteration. This can be done in different 
ways. The first and obvious one consists of the following 
steps: Solve the Kohn-Sham equations with an approx- 
imation to the OEP, e.g., do a KLI calculation. Use 
the resulting (p^ and v Kca to construct the right-hand 
side of Eq. l(2*T|) . Solve, e.g., by conjugate gradient 
iteration^S, for the orbital shifts. Insert the tjj* a into 
Eq. (|23|) to obtain an improved approximation for v xc<T . 
Then self-consistently solve the Kohn-Sham equations 
again for fixed u XCCT to obtain a new total energy and 
new eigenvalues. Repeat these steps until convergence is 
achieved. 

This simple procedure only involves the occupied or- 
bitals and their shifts. Instead of an integral equation, it 
only requires solving partial differential equations, which 



can easily be done. We thus gain all the advantages 
of avoiding the explicit evaluation of the response func- 
tion that are known from other areas of linear response 
theorji" 

In practice, a small modification of this algorithm is 
useful. After w XCCT has been updated, one need not di- 
rectly go back to the Kohn-Sham equations. Instead, one 
can solve Eq. I|21l) again, keeping the left-hand side and 
the ifiia- fixed but evaluating the right-hand side with the 
new v xcr7 . Using this additional (ip* a , w XC( j)-cycle speeds 
up the convergence, and Eq. (|21|l is easier to solve than 
Eq. Q since it is just a differential and not an eigenvalue 
equation. We discuss numerical aspects of the iteration 
and its convergence in detail in Appendix [CJ Here we 
just point out that the iteration converges very quickly: 
Starting , e.g. for the Be atom, from the KLI approxi- 
mation and using 5 cycles per iteration, the total energy 
is correct after the first iteration and all eigenvalues after 
5 iterations. 

The iterative approach is also accurate. One criterion 
that indicates an accurate OEP can be inferred directly 
from Eq. QJ. The function 



E 



CW^( r ) + c - c - 



(25) 



must vanish for the exact OEP. Therefore, the smaller 
the maximum value S max that S(r) takes in a given re- 
gion of space, the more accurate is the potential. For the 



KLI potential of the Be atom, we find S n 



0.03 a n 



After 5 iterations with 5 cycles, S max has fallen below 
10~ 6 (2q 3 . A second accuracy criterion is the exchange 
virial relation£i*Z£. It is slightly violated by the KLI ap- 
proximation with a relative error of about 1%. By the 
time our iteration has refined all eigenvalues to 0.0001 
hartree (0.1 mhar) accuracy, the error in the virial re- 
lation has been reduced to 10~ 4 %, and further iteration 
reduces it to 10 _7 %. The method is thus at least as ac- 
curate as the direct solution of the integral equation" 

However, applied to finite systems, Eq. i|23|) has an un- 
pleasant feature: The second term under the sum on the 
right-hand side must be divided by the density. While 
this is not a problem theoretically, it complicates the nu- 
merical evaluation of Eq. I|23[) for finite systems. Far- 
away from the system the density decays exponentially, 
and numerically dividing by a rapidly vanishing func- 
tion introduces inaccuracies. These inaccuracies in the 
asymptotic region can magnify during the iteration and, 
if not taken care off, can spoil the whole scheme. In 
Appendix [U] we discuss this problem in greater detail, 
together with strategies to overcome it. Here we want to 
focus on how it can be completely avoided. 

For many numerical problems, iterative solution strate- 
gies can be designed in different ways'. Therefore, there 
is no reason to believe that Eq. 1231) is the only way that 
the information contained in the orbital shifts can be 
used. In fact, Eq. suggests a very simple, prag- 

matic alternative: Since S a (r) is an indicator for the er- 
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ror inherent in a given approximation to the OEP, we 
improve the given approximation by adding the corre- 
sponding So- (r) to it: 

<c7(r) = <Or)+cMr). (26) 

The real parameter c > with dimension energy times 
volume is introduced because So-(r) is not an exact rep- 
resentation of the error but just an estimate. Since it is 
not guaranteed that the w"™(r) from Eq. i|26[l will obey 
Eq. Q24JI. we enforce it explicitly by adding the constant 
u xcN ^ - J v* Na (rK™(r)^ CT (r) d 3 r to v«™ each time 
it has been updated via Eq. I|26|) . (See Ref. |2J for an 
intuitive explanation of the same idea.) 

Eq. I)2(S|I . together with adding the constant, can re- 
place the update of w XC(T via Eq. J23J), and all the other- 
steps described in the three paragraphs following Eq. I|24|) 
stay the same. As shown below, introducing c can be un- 
derstood as replacing by a constant a term that requires 
dividing by the density. Doing so, we drop some informa- 
tion from the error feedback, and therefore the iteration 
converges somewhat slower than the direct iteration of 
Eq. (|2*5|) . But with 10 (^1*^, « XC(r )-cycles per iteration, 
the total energies of the atoms Be, Ne, Na, Mg, and 
Ar converge within 0.1 mhar in 2, 2, 2, 5, and 7 itera- 
tions, respectively (and a larger number of cycles further 
reduces the number of necessary iterations). The com- 
putational effort is therefore only moderate. The OEP 
energies from our solution are identical to the ones found 
from the considerably more involved direct solution of 
the integral equationS&iSi, and the iteration is as accu- 
rate (see Appendix ICl for tests). 

This approach can be motivated further by rewrit- 
ing Eq. I|21|) such that the left-hand side becomes 
(ft 2 / (2m)V 2 + £ ia )%l)* a and inserting the result into the 
right-hand side of Eq. I|22[) to obtain 

<w(r) = <w(r) + ^^(r). (27) 
2n CT (r) 

This is a trivial identity if </?j CT and tp* a are the orbitals 
and shifts corresponding to the true OEP. However, if 
(fiicr and have been calculated from an approximation 
to the OEP, then the second term on the right-hand side 
does not vanish and represents an error term. Its most 
important part is the function S a (r) which by definition 
contains the information about where the OEP equation 
is violated. By subtracting the error term from the ap- 
proximation, we get a better approximation. But we only 
want to subtract the part that is proportional to w XC(T , 

< c 7(r) = - -X°L) S< ,(t), (28) 

because for stability reasons we do not want to change 
Vxctr too drastically in each iteration. (Subtracting the 
full term using vks would, e.g., introduce a very different 
length scale and the nuclear singularity into the iteration. 
Besides, in the neutral systems of interest here, v xc dom- 
inates the asymptotics of uks-) We have tested Eq. I|28|) 



for spherical atoms, and it converges nearly as quickly 
as the direct iteration of Eq. 123|) . But since Eq. i|28|) 
also requires dividing by the density, the simpler itera- 
tion based on Eq. I|26l) is more practical, in particular for 
the cluster calculations. Eq. I|28|l provides an interpreta- 
tion for the constant c: It takes the role of an average 
value for w°J^(r)/(2n (T (r)). Since u xco - is negative and we 
chose c to be positive, the signs in Eqs. (12(i[) and l|28|) 
are opposite. In practice, the numerical value of c can 
be determined unproblematically by trial and error. For 
the systems we studied, we found c close to the values 
inferred from Eq. (|28|l : see Appendix IU1 



B. OEP for Na clusters 

Eq. I|26(l can very easily be employed in three- 
dimensional calculations and allows us to study the in- 
fluence that the exact exchange functional has on non- 
spherical systems. The electrons in atoms are rather 
strongly localized. Metal clusters, and in particular 
sodium clusters, have a very different electronic struc- 
ture with strongly delocalized valence electrons. They 
can thus be regarded as an "opposite" test case for the 
influence of exact exchange. Of particular interest is their 
static electric polarizability, which we will address in Sec- 
tion EH below. 

We first calculated the ground-state properties of small 
sodium clusters using the exact exchange functional to- 
gether with the exact OEP and with the KLI approxi- 
mation, and for comparison also with LDA for exchange 
only (xLDA) and LDA for exchange and correlation^. 
For reasons discussed in Section IVII we used the opti- 
mized cluster geometries and pseudopotential from Ref. 
F3l75l The results are summarized in Table H] For the 
smallest cluster Na2, which has one valence electron of 
each spin, the KLI approximation is exact. This is con- 
firmed by the exchange virial relation: Whereas for all 
other clusters the virial relation is slightly violated by 
the KLI approximation with a relative error of about 
0.5%, its relative error is only 10 -6 for Na2- The total 
energies of all larger clusters are lower with the OEP 
than with the KLI potential, as required. The abso- 
lute average difference of 0.5 mhar indicates that KLI 
is a rather good approximation for the total energies. It 
should also be noted, however, that the relative error is 
considerably larger for the clusters than for the atoms. 
The differences between the Kohn-Sham eigenvalues ob- 
tained from OEP and KLI are larger than the differences 
in total energy, and there is a clear pattern: As for the 
atoms, the OEP raises the lowest occupied eigenvalues 
and lowers the highest occupied ones. Thus, the energy 
range of the occupied Kohn-Sham spectrum is up to 6% 
narrower in OEP than in KLI. This might well influence 
time-dependent DFT calculations to which the eigenval- 
ues are an important ingredient. Comparing the results 
from the exact exchange functional to the xLDA results 
shows that, as expected, exact and local exchange lead 
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TABLE I: Absolute value of the total energy and lowest and highest occupied Kohn-Sham eigenvalue for small sodium clusters, 
in hartree, for exchange only OEP, KLI, and local exchange, and for local exchange and correlation^. For the two-valence- 
electron system Na2, the KLI potential is the exact OEP. 



Na 2 Na 4 Na 6 Na 8 Na i0 





-E — El = —£h 


-E — £l —£h 


-E — £l —£h 


-E — el —eh 


-E —£l —£h 


xOEP 


0.3768 0.1701 


0.7531 0.1762 0.1401 


1.1421 0.1861 0.1496 


1.5285 0.1954 0.1459 


1.9082 0.1956 0.1263 


xKLI 




0.7528 0.1776 0.1394 


1.1417 0.1871 0.1494 


1.5282 0.1963 0.1458 


1.9070 0.1965 0.1248 


xLDA 


0.3479 0.0874 


0.7047 0.1123 0.0717 


1.0810 0.1205 0.0824 


1.4619 0.1331 0.0826 


1.8257 0.1410 0.0664 


LDA 


0.4018 0.1135 


0.8166 0.1384 0.0976 


1.2511 0.1469 0.1085 


1.6926 0.1596 0.1085 


2.1169 0.1672 0.0926 



to considerably different eigenvalue spectra: Exact ex- 
change eigenvalues are more negative, and in particular 
the highest one is relatively much more negative. The to- 
tal energy in xLDA is less negative than that with exact 
exchange, reflecting the well known effect that the local 
approximation underestimates the magnitude of the ex- 
change energy. When local correlation is included, the 
local exchange error in the total energy is compensated 
to a good part. The eigenvalues, however, are still con- 
siderably smaller in absolute value even with correlation. 

Finally, we investigate the orbital shifts ip* a (r). This 
can be done in a transparent way for Na4: With two 
electrons of each spin, the KLI potential and OEP are 
already non-trivially different, yet we only need to look 
at two orbitals and shifts. The upper part of Fig.^shows 
the lower of the two Kohn-Sham orbitals as obtained from 
the KLI approximation and the OEP. The orbital is plot- 
ted for both potentials once along the y-axis and once 
along the z-axis. The ionic configuration of Na4 and our 
labeling of axes is shown in Fig. [21 Clearly, the orbital 
obtained within KLI is very similar to the orbital from 
the exact OEP. The second Kohn-Sham orbital which is 
not shown in Fig. ^ can be described as a p z orbital, i.e., 
it has a node in the x-y plane. The second KLI and OEP 
orbitals are also very similar. However, the conclusion 
that, therefore, the orbital shifts obtained from KLI and 
OEP should be very similar, too, is wrong. In the mid- 
dle part of Fig.Q]the first orbital shift is shown along the 
z-axis. The shift calculated from the OEP via Eq. 
shows a "dip" at the origin which is not seen in the shift 
corresponding to the KLI potential. And along the y-axis 
(lower part of Fig. , the shifts obtained from the KLI 
potential and the OEP are totally different: The KLI 
shift shows a pronounced central peak, whereas the OEP 
shift practically vanishes on the y-axis. The explanation 
for these somewhat surprising observations can be found 
directly in the OEP equation Eq. As mentioned be- 
fore, the x-y plane is a nodal plane for the higher of 
the two Kohn-Sham orbitals of Na4. The lower orbital, 
however, does not vanish in the x-y-plane. Eq. can 
therefore only be fulfilled in the x-y-plane if the orbital 
shift corresponding to the lower orbital vanishes in this 
plane. The KLI potential is only an approximation to 
the OEP, and therefore need not and, as seen in Fig. ^ 
does not fulfill Eq. Q in the x-y-plane. The iteratively 



constructed OEP, in contrast, makes the orbital shift van- 
ish in the plane, as it should. The lower part of Fig. 
thus once more confirms that we are indeed constructing 
the exact OEP. The reason for the pronounced difference 
between the OEP and KLI orbital shift is discussed in 
greater detail in Appendix IDl 

It remains to be explained why the orbital shift shown 
in the lower part of Fig. ^ vanishes perfectly between 
y=-7 ao and z=+7 ao, but still shows a tiny bump at 
larger distances, a remnant from the KLI orbital shift. 
The reason is our iteration based on Eq. \2ty ■ It corrects 
the potential quickly in the energetically-important re- 
gion where the orbitals and their shifts are appreciable. 
In the asymptotic region, however, fi a (j) and ip*a( r ) de- 
cay exponentially, and consequently the potential is mod- 
ified to a lesser extent. But why does this matter since 
it is believed that the KLI potential is asymptotically 
correct? The answer is that the KLI potential is asymp- 
totically correct everywhere except on nodal surfaces of 
the highest-occupied orbital that extend out to infinity. 
On such nodal surfaces, a small difference between KLI 
potential and OEP is observed even in the asymptotic 
region, and it is this difference which gives rise to the 
above mentioned "bump" . This aspect will be discussed 
in greater detail in the next section. Here we just briefly 
want to mention that, if needed, the iteration based 
on Eq. (|26[) can be modified to yield quickly the cor- 
rect asymptotic behavior also on the nodal surface: the 
constants v xc i a — u XC ia- which determine the asymptotic 
behavior (see below) only depend on the energetically- 
important region. Their values (and thereby the asymp- 
totic limits of the exchange potential) are thus accurately 
known after just a few iterations, even when the potential 
converges slowly in the asymptotic region. 



V. ASYMPTOTIC BEHAVIOR OF « x (r) 

The long-range asymptotic behavior of the exact ex- 
change OEP can be inferred from Eq. l|23l) . The first 
term in the sum over all occupied orbitals determines the 
asymptoticsi£i2&, since the highest-occupied orbital shift 
falls off faster than the highest-occupied orbital. The 
Kohn-Sham orbitals of a finite system decay exponen- 
tially like exp[— (— 2m£i cr /7i 2 ) 1 / 2 r], where r is the dis- 
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FIG. 1: The lowest Kohn-Sham orbital and corresponding 
orbital shift for Na4. Top: Full and long dashed: Orbital from 
OEP and KLI potential, respectively, along the y-axis. Short 
dashed and dotted: Orbital from OEP and KLI potential, 
respectively, along the z-axis. Middle: The corresponding 
orbital shifts plotted along the z-axis. Full: OEP, dashed: 
KLI. Bottom: Same orbital shifts along the y-axis. Full: OEP, 
dashed: KLI. See text for discussion. 



tance from the system's center. This is a consequence 
of the locality of the Kohn-Sham potential and means 
that each orbital's decay is dominated by its own orbital 
energy. Therefore, the sum asymptotically will be dom- 
inated by the slowest-decaying orbital, which is the one 
with the highest orbital energy. Since u X N a<T '—2? — e 2 /r, 





-0.05 - 


CD 


-0.1 - 


CD 










-0.15 - 






-0.2 - 




-0.25 - 




-0.3 - 



-30 




10 

yorz (a ) 



FIG. 2: The exact exchange potential of Na4 in hartree along 
the z-axis (full line) and the y-axis (dashed line) in bohr do- 
The potential goes to different asymptotic limits in different 
directions. 



v x (r) for a finite system goes to zero like — e 2 /r in all 
regions of space that are dominated by the highest- 
occupied orbital density nAr^(r) = |</Jjv CT ( r )| 2 - That the 
potential goes to zero and not some finite constant is due 
to Eq. J22J. However, there may also be asymptotic re- 
gions of space where n/v„ (r) vanishes, but the orbital den- 
sity of a lower occupied orbital M a does not. As pointed 
out in Ref. Hil this situation occurs whenever the high- 
est occupied orbital has a nodal surface that extends out 
to infinity. In these regions, one can go through all the 
arguments used to derive the asymptotic behavior and 
replace n^ a (v) by nj.f,(r). But the reference point for 
the potential was already fixed via Eq. (|24|) . We cannot 
choose a second one, therefore there is no equivalent of 
Eq. (|24|) for the M a orbital. Consequently, the potential 
will asymptotically go to C — v xc m„<x - u xcMa a, where 
in general C / 0. Thus, it must be concluded that the 
exact Kohn-Sham exchange potential can go to different 
asymptotic constants in different spatial directions. 

In Ref. |3(] this aspect of Eq. (|2"3*)l was first discussed, 
but its consequences for the Kohn-Sham potentials of 
real systems were not considered. In Ref. it was 
demonstrated that the effect indeed occurs for real sys- 
tems, but the numerical calculations were restricted to 
the "localized Hartree-Fock" approximation to the OEP. 
Finally, in Ref. 23 and in this paper, we evaluate the 
exact OEP including the asymptotic region for realis- 
tic, three-dimensional systems and confirm the surprising 
asymptotic behavior of the exact exchange potential. 

Fig. [21 shows the exchange-only OEP for the cluster 
Na 4 , the electronic structure of which was discussed in 
the previous section. The full line shows v x (r) along the 
z-axis. In this direction, the density is asymptotically 
dominated by the highest occupied orbital and the po- 
tential falls off like — e 2 /r, as expected. The dashed line 
shows v x (r) along the y-axis. Clearly, the potential falls 
of in the same way, but tends to a different constant 
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C. We find C = 0.0307 har. This non-vanishing value 
for C is of course found only in the x-y-plane. Mov- 
ing away from the origin at any low angle to this nodal 
plane, one eventually reaches a region where the highest- 
occupied orbital will dominate the density and v x (r) will 
tend to zero. A lower angle just means that the asymp- 
totic region starts further away from the center of the 
cluster. But nevertheless, since the exchange potential is 
a continuous function of r, the non- vanishing asymptotic 
constant has a substantial influence on i> x (r) even out- 
side of the nodal plane of the highest occupied orbital. 
To demonstrate this, we show in the lower half of Fig. 
13 the exact exchange potential for Na4, and the poten- 
tial from the local density approximation in the upper. 
For each point in the y-z-plane, the potential we cal- 
culated on a numerical grid is plotted as the "height". 
The LDA-potential mirrors the density closely, so the 
ionic cores which cause "bumps" in the valence electron 
density show up in the LDA potential. We also observe 
the well-known effect that the LDA potential decays very 
quickly: it vanishes on the boundary of the plotted re- 




21 

FIG. 3: The exchange potential "landscape" for Na4 in the 
y-z plane: w x (r) in hartree defines the "height" at each spa- 
tial point (y,z) for fixed x=0. Upper: exchange-only LDA 
potential; lower: exact exchange potential (OEP). Note the 
greater depth, slower decay, and pronounced "ridge" seen in 
the OEP. 



gion, whereas the exact exchange potential still takes a 
non-vanishing value there. But the most obvious differ- 
ence between the exact and the local exchange potential 
is that the former shows a pronounced "ridge" running 
along the y-axis. This ridge is a consequence of the non- 
vanishing asymptotic constant: The potential has to go 
up to reach its correct value on the nodal surface. We 
find qualitatively the same behavior for all the clusters we 
studied here. For Nag and Nas, there are two degenerate 
highest orbitals, and the OEP goes to C = 0.0378 and 
C = 0.0084, respectively, on the z-axis. (We choose our 
Cartesian coordinates such that the axes are the princi- 
pal axes of the cluster's tensor of inertia, and z is the axis 
whose principal value differs most from the average. All 
constants are in har.) The relatively small constant for 
the "magic" Nas reflects the fact that the highest orbital 
is nearly degenerate with the two lower-lying ones. For 
Naio, the situation is similar to N&i, with C = 0.0107 
in the x-y-plane. Finally, it is clear from Eq. I|23|) that 
the asymptotic constants are qualitatively correct in the 
KLI potential. But KLI and OEP orbitals differ slightly, 
as shown in the previous section. Consequently, we find 
constants that differ by about 0.001 har in all cases we 
investigated. 

The non-vanishing asymptotic constants are thus the 
rule, not the exception, and one may ask how this result 
can be true since there is a physical argument to explain 
the asymptotic behavior of the Kohn-Sham potential: An 
electron that wanders far out from a finite neutral sys- 
tem basically sees the attraction from the hole it leaves 
behind. Thus, the Kohn-Sham potential should fall of to 
— e 2 /r everywhere, without additional constants 76 . How- 
ever, appealing as this argument is, it cannot be applied 
to the Kohn-Sham potential. The Kohn-Sham potential 
is not a physical potential whose value we can measure. 
It is a mathematical construction, and therefore can have 
properties that physical potentials do not show. But of 
course the Kohn-Sham potential is a very clever math- 
ematical construction, designed to yield from a simple 
procedure an accurate description of complicated quan- 
tities we can measure, like the ground-state energy of a 
many body system. Therefore, it is important to un- 
derstand its features, in particular the surprising ones. 
The first and probably most extreme is the discontinu- 
ity of the Kohn-Sham potential as a function of electron 
number^S* 7 ^. The existence of non-vanishing asymp- 
totic constants is a second one. 



VI. STATIC ELECTRIC POLARIZABILITIES 
FROM THE EXACT EXCHANGE FUNCTIONAL 

To demonstrate the importance of incorporating ex- 
act exchange into functionals, we investigate the static 
electric dipole polarizability. There has been a long- 
standing controversy about the influence of the long- 
range asymptotics of the Kohn-Sham potential and the 
self-interaction error on the optical properties of metal 
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clustera 74 ! 80 ! 81 ! 82 . 83 ! 84 ! 85 ! 86 ! 87 ! 88 . The static electric dipole 
polarizability has been of particular interest: It was one 
of the earliest observables to give insight into the struc- 
ture of clusters^ 9 -, its measurement gives information on 
the effects that govern cluster growth^, and it may 
be used as an indicator for whether the bulk limit is 
reached 9 -!. The polarizability of small sodium clusters 
has been calculated beforo 74 i 80 i 83 i 85 i 87 i 88 . It was found 
that, while local and semi-local density functionals yield 
polarizabilities in reasonable qualitative agreement with 
experiment, they do not lead to quantitative agreement. 
However, it has also been demonstrated 9 ^ 3 ^ that a fi- 
nite cluster temperature greatly influences the static elec- 
tric polarizability, mainly through the effect of thermal 
expansion*^. The clusters' temperature in the experi- 
ments is presently not accurately known, but it is es- 
timated to be high enough to account perhaps for all of 
the observed differenced 4 * 9 ^ 3 .. Nevertheless, it is of great 
practical interest to obtain accurate theoretical values for 
the polarizability because the polarizability could then be 
used as a "thermometer" to estimate experimental clus- 
ter temperatures, as discussed in greater detail by Kronik 
and co-workers^— 

Our exact-exchange OEP calculations allow us to an- 
swer some of the questions related to the influence of 
self-interaction errors and long-range asymptotics that 
have been raised in the past. In the present study, our 
focus is on investigating these two effects by comparing 
LDA and exact-exchange results. In order to sort out 
the effects as clearly as possible, we want the valence 
electrons to move in exactly the same external potential 
for all E xc functionals. Thus, we use the same pseudopo- 
tential and cluster geometries for different E xc . For an 
accurate comparison with experimental data, a consis- 
tent exact-exchange pseudopotentia! 20 - 95,96 97 should of 
course be used and the cluster structures re-optimized. 
But the comparison with experiment would also require 
inclusion of correlation effects. Since a correlation func- 
tional generally compatible with full exact exchange is 
not yet available, these aspects will have to be addressed 
in future work. 

We calculated the mean static electric dipole polariz- 
ability a = tr(cc)/3 for the neutral sodium clusters with 
even electron numbers between 2 and 10, using the LDA- 
optimized geometries of Refs. 1740751 The polarizability 
tensor a was obtained from the electric dipole moment 
/ij(F) = —e J rjn(r,F)d 3 r by evaluating the definition 
aij = lim^^o d^ij/dFi, i,j = x,y , z, for a small but fi- 
nite electrical field F (see Ref. uM for details, e.g, on the 
role of the ionic cores) . 

Table [H] summarizes the results of our calculations. 
We find that, for exact exchange and no correlation, the 
polarizabilities calculated using the exact OEP are con- 
sistently lower than the ones calculated from the KLI ap- 
proximation. This is in agreement with our expectations, 
since the OEP yields lower total energies, i.e, binds the 
valence electrons more strongly than the KLI potential, 
and the valence electron density is thus harder to displace 



by the electric field. It should be noted, however, that the 
difference between OEP and KLI polarizabilities is on the 
order of 1 %. By repeating our calculations for different 
electric field strengths and on different Cartesian grids 
with spacings between 0.5eto and 0.8eto and between 65 
and 129 point in each spatial direction, we estimated the 
numerical accuracy of our calculations to be somewhat 
better than 1%. Thus, the inaccuracies that are intro- 
duced by using the KLI approximation are close to the 
numerical limits, showing once again that, for the present 
systems and observables, KLI is a good approximation to 
the OEP. 

Compared to exchange-only LDA, the OEP yields po- 
larizabilities which are on average a noticeable 7% lower. 
Again, this result can readily be understood: Local ex- 
change makes a large self-interaction error, i.e., overes- 
timates the electron-electron repulsion. This leads to a 
density that is too delocalized, which in turn results in 
too high a polarizability. The self-interaction error de- 
creases with increasing electron delocalization, and the 
valence electrons in sodium are nearly free electrons. 
Therefore, we expect the largest self-interaction errors 
and differences between local and exact- exchange polar- 
izabilities for the smallest clusters. This is confirmed by 
Table ILT1 However, the fact that the difference falls from 
Na2 to Nas, but increases somewhat from Na§ to Naio, 
shows that this argument seems to hold for the over- 
all trend, but details in the electronic configuration also 
seem to matter. 

Finally, comparing exact-exchange OEP to LDA for 
exchange and correlation, we find the surprising result 
that the polarizabilities are nearly the same. How can 
this be understood? The self-interaction error for local 
exchange and correlation together is smaller than for lo- 
cal exchange alone. So if inclusion of local correlation 
for the sodium clusters would lead to perfect cancella- 
tion of the self-interaction errors, then the fact that the 
polarizabilities are practically equal would not come as a 
surprise. However, the local correlation functional does 
not merely reduce the self-interaction error. Correla- 
tion functionals quite generally also describe short-range 



TABLE II: Mean static electric dipole polarizability a of Na 
clusters in A 3 for different approximations to E xc and v xc . 
Columns from left to right: exact exchange with exact OEP; 
exact exchange with potential in KLI approximation; LDA for 
exchange only; LDA—- for exchange and correlation; relative 
difference in % between local and exact-exchange polarizabil- 
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Coulomb effects that lead to stronger binding. Thus, we 
expect that adding an appropriate correlation functional 
to the exact exchange will bind the valence electrons more 
strongly, yield a more compact density and thus lower 
the polarizability. Qualitatively, our results thus indi- 
cate that the LDA overestimates the polarizability. This 
conclusion, of course, is restricted to our special situation 
of fixed external potential. In future work we plan to in- 
vestigate this question further by combining a newly de- 
veloped, self-interaction-free mcta-GGA for correlation 98 
with the exact exchange functional. 



VII. SUMMARY AND CONCLUSION 

In summary, we have presented a new proof for 
the OEP equation that highlights the relationship be- 
tween different density functionals (e.g., Hartree-Fock 
and exchange-only OEP) that can be defined from a given 
orbital expression for the energy. The proof explains why 
the complicated OEP equation can be cast into the sim- 
ple form of a vanishing density shift. We discussed how 
the exact OEP can easily be calculated from only the oc- 
cupied Kohn-Sham orbitals and their first-order shifts. 
The exact exchange potential was calculated for (all- 
electron) closed-shell atoms and three-dimensional (pseu- 
dopotential) sodium clusters. For the first time, not only 
the total energy but also Kohn-Sham eigenvalues were 
accurately calculated for three-dimensional systems from 
the exact OEP. We discussed different indicators for the 
accuracy of an OEP solution. The KLI approach is a 
rather good approximation for the ground-state proper- 
ties of atoms and small sodium clusters. The relative 
errors it introduces, however, are considerably larger for 
the clusters than for the atoms. We further investigated 
the asymptotic behavior of the exact exchange poten- 
tial of Kohn-Sham theory. For the non-spherical clus- 
ters, it shows pronounced "ridges" which are missing in 
the potentials of widely-used approximations like LDA. 
Finally, we calculated the static electric dipole polariz- 
ability for several sodium clusters and compared exact- 
exchange OEP to exact-exchange KLI and to LDA. The 
advantages of self-interaction-free exact exchange were 
demonstrated, and the need for a correlation functional 
to go with exact exchange stressed. 

Our method greatly facilitates the self-consistent use 
of orbital functionals like the exact exchange energy in 
Kohn-Sham calculations. Orbital functionals can now be 
used fully self-consistently with moderate effort, and the 
accuracy of approximations to the OEP, such as KLI, 
can be checked. Due to the correct long-range asymp- 
totics, functionals that use full exact exchange appear 
very promising for time-dependent DFT. The great influ- 
ence that exact exchange has on the Kohn-Sham eigen- 
values, which play a prominent role in time-dependent 
calculations, was demonstrated for the metal clusters. 
The locality of the Kohn-Sham potential is particularly 
beneficial in this context and makes Kohn-Sham theory 



superior— to theories with orbital-dependent potentials, 
e.g., Hartree-Fock. Therefore, we believe that orbital 
functionals and the OEP method will play a prominent 
role in future functional development. 
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APPENDIX A: THE ORBITAL ENERGY 
DENSITY FUNCTIONAL E° Th ' a [n] 

Our particular definition of the orbital functional 
E° lh ' a [n] was chosen to avoid problems that other rather 
natural-looking definitions may have. An intuitively 
plausible definition, e.g, could have required the orbitals 
to be orthogonal. However, it is known that certain or- 
bital functionals, e.g., the SIC functional, are mini- 
mized by non-orthogonal orbitals. Requiring orthogonal- 
ity in the definition of E° rh ' a [n] for the SIC would require 
non-trivial modifications of Eq. I|14(l . e.g., introducing 
off-diagonal Lagrangian multipliers. On the other hand, 
simply dropping the orthogonality condition and requir- 
ing the orbitals to be just the ones minimizing the energy 
could lead to a "bosonic" solution, with all electrons oc- 
cupying the lowest orbital. Therefore, in the definition of 
E° lh ' a [n] the orbitals are required to be linearly indepen- 
dent. They are further required to be extremizing ones 
and not just the ones yielding the lowest possible energy 
for a given density, because the latter criterion could be 
fulfilled by a "nearly bosonic" solution, with all electrons 
occupying orbitals that are extremely similar to the low- 
est one but with a tiny contribution added to each one 
of them to make them linearly independent. Our defini- 
tion of E° rh ' a [n] raises a question similar to the questions 
of w-representability which are still a subject of current 
research in Kohn-Sham theory— In our case the ques- 
tion is: Can we find a set of orbitals that extremizes Eq. 
@ for any given density? However, we believe that for 
our present purposes this question does not pose a prob- 
lem because we are interested in the minimizing density, 
which we know can be obtained from Eq. 114|l . 

Eq. (|14|) starts from a given external potential v(r) 
and finds the orbitals and density belonging to it. If 
instead we want the (/9° rb,Q! [n] corresponding to a given 
density n(r) , we must add to v(r) in Eq. I|14|) a Lagrange- 
multiplier Ai>([n]; r) which enforces the constraint of Eq. 
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APPENDIX B: PROOF THAT THE MINIMIZING 
DENSITIES FOR E° Th ' a [n] AND E° EP ' a [n] AGREE 
TO LINEAR ORDER IN a 

In Eqs. I|12|) an d (|13|l . we denned the functional 
E° EP ' a {n] and E° lh < a {n] by starting from Eq. © and 
its discussion. Assuming that -E° rb ' a [n] is analytic in a, 
we consider the power series expansion 



E° lh ' a [n] = E° EP ' a [n] + aA 1 [n] + a 2 A 2 [n] + .. 
The analog of this expansion for the densities is 

n orb ' Q (r) = n OEP ' a (r) + cmi(r) + a 2 n 2 (r) + ... 
We also argued in Section ITT1 that 
Ai[n] = 0. 

Now rewrite the Euler equation for n orb ' Q (r): 
5E°* h ' a h 



6n(r) 

SE? EP ' a \n] 







Sn(r) 
6E° EP ' a \ 



i SA 2 [n] 



Sn(r) 
SE? EP ' a \n] 



Sn(r) 

+ 0{a 2 ) 



(Bl) 
(B2) 
(B3) 

(B4) 
(B5) 
(B6) 



i OEP -°+an 1 +. 



Sn(r) 



n OBP,o 



6 2 E OEP„ 



n 



Sn(r)Sn(r') 



ni(r')dV + C>(a 2 ) (B7) 



+a 



Eq. (|Bg) follows from Eqs. (JHTJ and $5J$. In Eq. JEM)) we 

use Eq. l)B2j) . and we can focus on the terms linear in a 
because each order of a must vanish separately. Eq. I|B7|) 
finally follows from functional Taylor expansion. The 
first term in (|B7|) is just the Euler equation for n OEP,a (r) 
and therefore vanishes. So 



5 2 E 



OEP,q[ 



5n(r)6n(r') 



ni(r')dV = 0, 



(B8) 



and Eq. (|B8|) must be fulfilled for all a. Since the double 
functional derivative in Eq. (|B8|I depends upon a, while 
ni(r') does not, Eq. I|B8(I can only be fulfilled by 



ni(r) = 0, 



(B9) 



i.e, the minimizing densities n orb,Q (r) and n OEP,Q (r) 
agree to linear order in a. 

We can compute n.i(r) from first-order perturbation 
theory. Since, from Eq. ljB2|) . 



cmi(r) = n orb < Q (r) - n OEP ' Q (r) + 0{a 2 ), (BIO) 

ni (r) can be computed by evaluating the density change 
that is associated with going over from the OEP- 
functional to the orbital density functional and neglecting 



all orders in a beyond the linear one. This density change 
can be calculated from the corresponding change in the 
single-particle orbitals. To this end, we rewrite Eq. Ijl4|l 
with the help of Eqs. ljT7J) and ((THJl in the form 



V 2 + v(r) + v H (r) + av xc (r) 
2m / 

i r / ^ / \ i orb, a / \ orb, a orb, a / \ /t-> i -i \ 

+a[u xci (r) -v xc (v)]ip l (r)=e i ip, (r).(Bll) 



Clearly, if we neglect the term in brackets in the second 
line of Eq. (|B11|) , we are looking at the Kohn-Sham sys- 
tem. Thus, we can calculate the shift that the orbitals 
undergo when we go over from the Kohn-Sham to the 
orbital density functional by treating a[u XC i(r) — v xc (r)] 
as a perturbation on the Kohn-Sham system. In non- 
degenerate first-order perturbation theory (see SectionHTI 
for remarks on the degenerate case), the shifted orbital 
is 



pert, a/ \ _ 



(r) = ^(r) + 



(<pf\a [u xcl - v xc ] \<pf) 



3 = 1 



Si 



where yf(r) and Si are the Kohn-Sham orbitals and 
eigenvalues corresponding to the minimizing density 
n p,Q (r). The important fact to be noted here is that 
the perturbation is linear in a. Thus, higher orders of 
perturbation theory will also be of higher order in a. 
Since we want to calculate the change in density only to 
first order in a, only the orbital shift in first-order pertur- 
bation theory needs to be taken into account. We insert 
Eq. (HUll into Eq. ijBTOl) . 



N 



N 



onx(r) = ^^ crt ' a (r)^ pCrt ' a *W-E^( r )^*(r) 
i=l i=l 
JV 

= aJ2tf(r)<Ph*(r)+c.c. + 0(a 2 ), (B13) 



where ipf 1 (r) is defined as the second term on the 
right-hand side of Eq. (|B 1 2|l divided by a. Obviously, 
yfr 1 (r) = -ipi(r), and thus Eq. (|B9|l is exactly Eq. JTJ. 



APPENDIX C: NUMERICAL ASPECTS OF THE 
ITERATIVE CONSTRUCTION 

In this Appendix we summarize some of the more tech- 
nical details of the iterative solution of the OEP equation. 

We first discuss the convergence of the iteration based 
on Eq. (|23|l for the exact-exchange-only OEP. The test 
system here is the Be atom. The Kohn-Sham equations 
were solved on a logarithmic radial grid. Table [ffl] illus- 
trates how many iterations going back and forth between 
Eq. J3J and Eq. (|23|) are necessary for a given number of 
(ipia-i w XC cr)-cycles to converge the total energy and eigen- 
values to the OEP values with an accuracy of 0.1 mhar. 
The starting point is the KLI potential. The total energy 
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is correct after one iteration in all cases (the difference is 
small to begin with), but also the eigenvalues which are 
noticeably different in KLI and OEP converge quickly: 
With just 5 (ipicn ^xcoO-cycles per iteration, it takes only 
5 iterations to reach the exact OEP. We also note that for 
the second and later iteration steps, the self-consistent so- 
lution of the Kohn-Sham equations proceeds much faster 
than the first self-consistent solution that is needed to get 
the starting approximation. This has two reasons. First, 
we need to recalculate the u^da only once at the begin- 
ning of each new (V>| CT ) Vxca)- C y c ^ and not at each step in 
the self-consistent solution of the Kohn-Sham equations 
(since v xc is kept fixed during the latter). Second, the or- 
bitals and shifts from the previous iteration are already 
a good starting guess for the next one. Therefore, the 
whole scheme converges very rapidly. 

However, Eq. Q23JI requires dividing by the density, and 
this is a disadvantage if it is to be used for finite sys- 
tems. The first term under the sum is unproblcmatic. 
It is asymptotically dominated by the highest occupied 
orbital which also dominates the density. Thus, even if 
the density in the denominator has decayed to the extent 
that its numerical representation becomes highly inaccu- 
rate, the same inaccuracies show up in the numerator, 
and the correct result is still obtained. This cancella- 
tion does not occur in the second or gradient term, and 
the resulting inaccuracies can make the scheme unsta- 
ble. For spherical systems, the problem can easily be 
solved. In the asymptotic region, the KLI potential is 
already correct, therefore the second term need not be 
evaluated numerically there. And the asymptotic region 
for a spherical system can easily be identified by going 
out from the system's center and monitoring when the 
density starts to be dominated by the highest occupied 
orbital. In this way, calculations for spherical atoms can 
be done simply. 

But the accuracy of the method depends rather sensi- 
tively on the cut-off point: With a cut-off too close to the 
system's center, too much of the second term is lost; with 
a cut-off that is too far away, the iteration becomes unsta- 
ble. In the spherical, effectively one-dimensional calcula- 
tions, the numerical parameters can be chosen such that 
finding an appropriate cut-off radius is possible. Our real 



TABLE III: Number of iterations necessary to converge all 
eigenvalues of the Be atom to their exact OEP values E = 
-14.5724, ei = -4.1257, e 2 = -0.3092, for a given number of 
(ip* a , VxcoO-cycles. The total energy is correct after one itera- 
tion. Starting point is the KLI approximation (E = —14.5723, 
ei = —4.1668, £2 = —0.3089, all energies in hartree). Note 
that convergence is achieved with very moderate effort. 



(V , * t T>"xc,r)-cycles 


iterations 


(VC«xc CT )-cycles 


iterations 


1 


24 


5 


5 


2 


10 


10 


3 


3 


9 


20 


3 



interest, however, is in non-spherical systems of possi- 
bly complex shape. We found that for three-dimensional 
(3D) systems, the technical difficulties are more severe. 
The coarser grids that computational necessity forces on 
us lead to larger inaccuracies in the evaluation of the 
derivatives that appear in the second term and allow for 
less fine tuning of the cut-off point. Furthermore, the 
asymptotic region will start at different distances in dif- 
ferent directions and must be mapped out more carefully. 

As a first test case for a 3D implementation of our al- 
gorithm, we chose the jellium cluster Nas. It is a good 
test case because it is spherical, i.e., the OEP can be 
obtained from the established direct solution of the inte- 
gral equation for comparison^. But unlike the diverging 
atomic nuclear potential, jellium Nas can a l so be repre- 
sented on the 3D Cartesian grid which we use for a re- 
alistic, pseudopotential-based description of metal clus- 
ters. The Seitz radius {r s = [3/(47m)] 1 / 3 ) of the jel- 
lium background density was r s = 3.93 ao, correspond- 
ing to a sphere of radius 7.86 ao- We first calculated 
the KLI approximation and compared a spherical, one- 
dimensional calculation on a logarithmic grid to a 3D 
calculation (Cartesian grid, 65 points in each direction, 
spacing 0.7<xo, damped gradient iteration with finite- 
difference multigrid relaxation^) that does not exploit 
the spherical symmetry. Both yield identical energies and 
eigenvalues to within 0.1 mhar. This confirms the accu- 
racy of our 3D code. We then constructed the OEP itera- 
tively on the same grid for different cut-off radii between 
4 ao and 19 ao. Only with cut-offs close to 6 ao was the it- 
eration unproblematic and converged to the correct total 
energy E=-0.3735 har. However, the lowest eigenvalue 
was still 0.4 mhar above its correct value 0.2081 har, i.e., 
halfway between the OEP and KLI values. Obtaining 
the OEP via Eq. I|23|) for the real, non-spherical clusters 
turned out to be even more complicated. The situation 
can be improved with finer grids, but using Eq. I)23|l in 
3D calculations is rather tedious. 

The iteration based on Eq. (|26|) does not suffer from 
these problems. The convergence is slower than for the 
first method, but the computational effort is still very 
moderate: For the test system jellium Nas, the total en- 
ergy and all eigenvalues converge within 5 iterations if 
10 cycles are used. With Eq. J2EJ), the 3D calculation 
now also yields exactly the same energy and eigenval- 
ues as the spherical, direct solution of the OEP integral 
equations^. During the iterations, the relative error in 
the exchange virial relation falls by two orders of mag- 
nitude compared to the KLI approximation, and S'max 
decreases to 2 x 10~ 6 . Further iterations reduce both 
quantities further. Thus, the simple iteration very accu- 
rately yields the exact OEP also for 3D systems. 

The parameter c used in the iteration was chosen by 
trial and error. If c is chosen too large, the iteration 
diverges, and with c too small, it does not converge as 
quickly as possible. In all cases we studied, finding ap- 
propriate values for c was unproblematic and there was 
a range of suitable values. We used c's between 3 har 
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and 0.75 har Oq for the atoms, and c—30 har ajj for all the 
cluster calculations. Thus, c was similar for similar sys- 
tems, i.e., c's of about 1 were appropriate for the atoms 
with their singular nuclear potential, and larger c's could 
be used for the non-singular pseudopotential calculations. 

These values can be understood on the basis of Eq. 
(|28p. For Na4, e.g., — v XCT (r)/(2n cr (r)) takes values be- 
tween 30 har Oq and 40 har a§ in the interior region of 
the cluster. A yet simpler and charming argument can 
be made based on the bulk limit. For a homogenous sys- 
tem, —Vx/n = (97r/4)( 1 / 3 Vs har ao » 1.92 harao- For 
bulk sodium, r s = 3.93 ao, so —v x /n w 30 haraj], and for 
atoms, the average density is about r s — lao, so —v K /n 
is about 2 hara^. These estimates agree nicely with the 
values we found by trial-and-error. Finally, we want to 
remark that the iteration could also be used with some 
suitably chosen c(r) , but for our present purposes a con- 
stant was satisfactory. 

APPENDIX D: ON THE DIFFERENCE 
BETWEEN KLI AND OEP ORBITAL SHIFTS 

Looking at Fig.[TJ one might ask the following question: 
The KLI and OEP orbitals are rather similar. Therefore, 



in Eq. i|21|) the operator on the left-hand side and all 
quantities on the right-hand side should be similar, too. 
But then, also the orbital shifts that are calculated from 
Eq. (|21l) should be similar. So how can it be that we 
observe the pronounced difference between the KLI and 
OEP orbital shift ipi (r) shown in the lowest part of Fig. 
Iff 

The answer to this equestion is that although the 
Kohn-Sham orbitals obtained from the KLI potential and 
the OEP are similar, the respective right-hand sides in 
Eq. I|21l) nevertheless show a qualitative difference. The 
higher Kohn-Sham orbital <^2(r) has a nodal surface in 
the x-y plane, i.e., its orbital density vanishes there. 
Therefore, in this plane the KLI potential is just given 
by Wxi(r) — (Sxi — u x i). If this is inserted for v x (r) in the 
right-hand side of Eq. Q21[l. the right-hand side vanishes. 
But for the OEP, there is an additional contribution to 
i> x (r) from the term of Eq. 1)23(1 that involves the deriva- 
tives of the Kohn-Sham orbitals and their orbital shifts, 
and this term makes a non- vanishing contribution in the 
x-y plane. The qualitative difference between OEP and 
KLI orbital shift is thus a direct consequence of a quali- 
tative difference of the respective right-hand sides in Eq. 
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